rm(list=ls())
setwd("E:/data/cooperation/finance_research/innovation/data")
library(foreign)
library(xlsxjars)
library(xlsx)
library(readstata13)
library(stringr)
library(ggplot2)
library(reshape2)
library(gridExtra)
library(randomForest)
library(e1071)
library(party)
library(mapchina)
library(sysfonts)
library(showtextdb)
library(showtext)
library(tidyverse)
library(rlang)
library(ggplot2)
library(sf)
library(dplyr)
library(RColorBrewer)
library(showtext)
##
da<-read.dta13("co2_finance.dta")


a<- floor((da$year-1999)/10)
a<-as.factor(a)
da<-data.frame(da,year1=a)
q1<-ggplot(da,aes(x=fscale,y=per_co2))+
  geom_point()+
  geom_smooth(method="lm",formula = y ~log(x)+x^3,size=1,color="black")+
  theme_bw()+
  #xlim(1,100)+ 
   ylim(0,20)+
  #scale_color_gradient(low = "gray",high = "black")+
  #scale_fill_manual(name="Year",values=c("gray40", "red"),labels=c('2000-2008','2009-2017'))+
  theme(title=element_text(size=12),text=element_text(size=12,family="serif"))+
  labs(y="Carbon dioxide emissions per capita",x="Financial development",title="(a) Scatter plot and fitted curve")+
  theme(panel.grid=element_blank(),panel.border=element_blank(),plot.title = element_text(hjust = 0.5),axis.line = element_line())
##
da$per_co2=da$per_co2/10
a<- floor((da$year-1999)/10)
a<-as.factor(a)
da<-data.frame(da,year1=a)
ggplot(da,aes(x=fscale,y=per_co2))+
  geom_point(aes(shape=year1,colour=factor(year1),fill=factor(year1)),size=1.2)+
  geom_smooth(method="lm",formula = y ~log(x)+x^3,size=1,color="black")+
  theme_bw()+
  #xlim(0.1,1)+ 
  ylim(0,2)+
  scale_color_manual(name="Year",values=c("grey30", "grey2"),labels=c('2000-2009','2010-2017'))+
  scale_fill_manual(name="Year",values=c("grey30", "grey2"),labels=c('2000-2009','2010-2017'))+
  scale_shape_manual(name="Year",values=c(1, 17),labels=c('2000-2009','2010-2017'))+
  #scale_linetype_manual(name="Year",values=c(4, 1),labels=c('2000-2009','2010-2017'))+
  theme(title=element_text(size=10),text=element_text(size=10,family="serif"))+
  labs(y="Carbon dioxide emissions per capita",x="Financial depth")+
  theme(panel.grid=element_blank(),panel.border=element_blank(),plot.title = element_text(hjust = 0.5),axis.line = element_line())



#
da$dp=round(da$fscale)
dp<-1:8
da_1=data.frame(year=dp)
da_1$mean<-0
for (i in 1:8)
{
  da_s<-da[which(da$dp==(i)),]
  da_1$mean[i]<-mean(da_s$per_co2,na.rm = T)
}



q2<-ggplot(da)+
  geom_boxplot(aes(x=dp,y=per_co2,group=dp),color="black")+
  geom_point(aes(x=dp,y=per_co2),alpha=0.3,position = "jitter",size=0.8)+
  geom_point(data=da_1,aes(x=dp,y=mean),shape=17,size=3)+
  geom_line(data=da_1,aes(x=dp,y=mean),size=1.5,linetype="dotted")+
  theme_bw()+
  #ylim(1,20)+
  scale_x_continuous(breaks=seq(1,8,1))+
  labs(x="Financial development",y=NULL,title="(b) Boxplot and mean curve")+
  theme(title=element_text(size=12),text=element_text(size=12,family="serif"))+
  theme(panel.grid=element_blank(),panel.border=element_blank(),plot.title = element_text(hjust = 0.5),axis.line = element_line())


grid.arrange(q1,q2,ncol=2)

##
df3 <- china 

df3 <- df3 %>%
  group_by(Name_Province) %>%
  summarise(geometry = st_union(geometry))
a<-df3$Name_Province
#write.csv(a,file="中国各省.csv")
b<-read.csv("energy_tech.csv",header=T)
df3$Fossild<- b$energy_d
ggplot(data = df3) +
  geom_sf(aes(fill = Fossild)) +
  scale_fill_distiller(palette = "RdGy",name="Regional distribution of fossil energy dependence") +
  theme_bw() +
  labs(title="(b) Area division")+
  theme(text=element_text(size=10,family="serif"))+
  theme(legend.position = "bottom",plot.title = element_text(hjust = 0.5,size=9))



##
df3 <- china 

df3 <- df3 %>%
  group_by(Name_Province) %>%
  summarise(geometry = st_union(geometry))
a<-df3$Name_Province


b<-read.csv("energy_tech.csv",header=T)

df3$Fossild <- b$tech_d
ggplot(data = df3) +
  geom_sf(aes(fill = Fossild)) +
  scale_fill_distiller(palette = "RdGy",name="Technology innovation (time mean)") +
  theme_bw() +
  labs(title="(b) Area division")+
  theme(text=element_text(size=10,family="serif"))+
  theme(legend.position = "bottom",plot.title = element_text(hjust = 0.5,size=9))


##中国地图

setwd(mywd)
# 载入所需包
library(sf)
library(dplyr)
library(ggplot2)
# 读入绘图需要的地图数据:省界、国界、九段线
provinces
boundary
lines9
# 以上地图没有定义坐标系统,给它们统一配上beijing54地理坐标系统，其EPSG编码为4214
st_crs(provinces)
st_crs(boundary)
st_crs(lines9)
# 读入属性数据(这里为各省份2010年人均GDP),并与provinces对象进行合并
pcgdp
provinces
# 为实现分层效果,将2010年人均GDP数值切段为因子
provinces$cut2010

ggplot() + geom_sf(data=provinces,  mapping=aes(fill=cut2010)) +
  
  geom_sf(data=boundary) +
  
  geom_sf(data=lines9, size=1.2,  colour="grey30") +
  
  coord_sf(crs=21458) +
  
  theme(legend.position=c(0,0),
        
        legend.justification=c(0,0),
        
        legend.background=element_blank()) +
  
  guides(fill=guide_legend(title="人均GDP2010"))



x=seq(0,9,0.01)
y=3.2558*x-0.3868*x^2-0.5

b=-0.3868
c=3.2558
a<-4.56
a<-(-c)/(b*2)
3.2558*a-0.3868*a^2-0.5

da<-data.frame(x,y)
da$color<-0
da$color[which(da$x>=4.208635)]<-1

ggplot(da,aes(x=x,y=y))+
  geom_line(size=0.9)+
  theme_bw()+
  theme(title=element_text(size=10),text=element_text(size=10,family="serif"))+
  geom_hline(yintercept =6.351237,lty="dashed")+
  annotate(geom = "segment",x=4.208635,
           xend = 4.208635,y=6.351237,yend = -Inf,
           lty="dashed")+
  geom_point(aes(x=4.208635,y=6.351237),size=3,shape=17,
             color="black",alpha=0.9)+
  annotate(geom = "text",x=4.438635,y=1.8,
           label=round(4.208635,2),size=3,
           vjust=6,color="black")+
  coord_cartesian(clip = "off")+
  xlim(1,8)+
  ylim(1,8)+
  labs(y="Carbon dioxide emissions per capita",x="Financial development")+
    #scale_x_continuous(breaks=seq(1,9,2))+
    #scale_y_continuous(breaks=seq(0,8,2))+
  theme(panel.grid=element_blank(),panel.border=element_blank(),plot.title = element_text(hjust = 0.5),axis.line = element_line())




